#!/usr/bin/perl
use strict;
use warnings;
use File::Copy;
use File::Basename;

################################ Overview ########################################

#Processes yseq subdirectories of one or more YCR directories of a cross
#combining the raw haplotype calls, carrying out QC for markers and strains across 
#the whole dataset and inferring missing values using multiple approaches.




#Will infer the genotypes of dead spores coming from tetrads the other members of
#which are sequenced


##################################################################################

#Written by Gareth A. Cromie, 1st November 2011

################################## Usage #########################################

# perl whole_cross_haplotype_call_V1_4.pl <cross_directory><cross><strain_file><ycr_dirs>

#<cross_directory>=directory holding full cross haplotypes
#<cross>= name of strains used in cross in format "YO490_YO502" will create subdirectory
#	  of cross_directory for the relevant cross (if not already existing
#	  each time script run for a given cross will create a new "cross_haplotype_X"
#	  subdirectory of the cross subdirectory
#<strain_file>=LIMS YPG csv output file allowing dead spores associated with
#		tetrads having sequenced members to be associated with the tetrad and
#		have their genotype inferred (allowing linkage analysis for phenotype
#		"DEAD")
#<ycr_dirs>=the inividual YCR directories holding the YSEQ raw haplotype files



##################################################################################
##################################################################################
################################### Main #########################################
##################################################################################
##################################################################################


##################################################################################
######################### Initialize Directory Structure #########################
##################################################################################


my $version="full_cross_haplotype_call_V1_4.pl";
my $usage="$version <cross_dir><cross><strain_file><ycr_dirs>";
unless (@ARGV>1){die "perl $usage\n";}
my ($cross_dir, $cross, $dead_file, @ycr_dirs)=@ARGV;
my $r_path="/tools/R-2.13.1/bin/R";



#Create an output directory structure
unless (-d "$cross_dir/$cross"){mkdir "$cross_dir/$cross" or die "Could not create haplotype master directory:\n$!\n";}
 
my $count=1;	#Use count to create unique combined haplotype directory name
while (-d "$cross_dir/$cross/cross_haplotype_$count"){++$count;}
my $outdir="$cross_dir/$cross/cross_haplotype_$count";
mkdir $outdir, 0755 or die "Could not create current cross haplotype output directory $outdir:\n$!\n";

my $label="${cross}_cross_haplotype_${count}";	#Label identifying cross-iteration


#Create logfile
open (LOG, ">${outdir}/log.tsv") or die "Could not open log file ${outdir}/log.tsv\n";
print LOG "perl $version\t@ARGV\n";
close LOG;

##################################################################################
################# Create Hashes by Strain, Marker and Tetrad #####################
##################################################################################

#Holds information about each marker instance: strain, chromosome and position as 
# keys to array of tag, plate, tetrad, spore, p1 confidence, p2 confidence, total
#informative reads, allele call,P1 perfect match reads, P2 perfect match reads, 
#perfect match allele call, yseq number, cross number.
my %markers=();	 

#Holds information about tetrads with plate, tetrad and strain as keys to value 1
#if tetrad has consistent marker segregation (good) or 0 if not
my %tetrads=();	

#Holds information about all markers ever called:chromosome, position as keys to 
#marker assessment hash: word "type" as keys to BA (bi-allelic- parent specific 
#reads expected from both parents) or MA (mono-allelic, restriction site present 
#in only one parent) and allele (0,1,2,3 etc) as key to counts of that allele
#across all strains (filled in by assess_markers subroutine)
my %positions=();

#Keep track of strain information:strain as key to array of tag, plate, tetrad,
#spore, YSEQ and YCR (cross)
my %strain_info=();	
			

##################################################################################
############### Combine Raw Haplotypes From all Strain-YSEQs #####################
##################################################################################

#Get yseq directories for all crosses and combine all of the raw haplotype file 
#data in the appropriate hashes
my @folders=();
foreach my $ycr (@ycr_dirs){
	@folders=(@folders, glob "$ycr/yseq/YSEQ*");
}

&combine_yseqs(\%markers, \%tetrads,\%positions, \%strain_info, \@folders);
		
##################################################################################
############################ Assess Strain Metrics ###############################
##################################################################################

#Assess each strain for quality metrics:strain as key to array:
#Total reads, % informative markers called as allele3 (likely mixed sample), 
#Hash holding ploidy for each chromosome (not 2 micron or mitochondria),
#Percentage 3:0 calls of P1 or P2  among the other 3 members of the same tetrad,
#	where all calls P1 or P2 (single low value when rest of tetrad high 
#	indicates single spore incorrectly assigned to tetrad)
my %assess=();
&assess_strains(\%assess, \%tetrads, \%markers, \%positions);

#Assess number of markers in each strain called  and number called as P1 or P2
&p1_p2_support(\%markers, "$outdir/initial_p1_p2_support.tsv");

#Output assessment results
open (ASSESS, ">$outdir/strain_assess.tsv");
print ASSESS "#Strain\tPlate\tTetrad\tSpore\tTotal_Reads\t%_Het\t%_Non_mendelian_seg_rest_tetrad";
print ASSESS "\t$_" foreach (sort keys %positions);
print ASSESS "\n";
foreach my $plate (sort keys %tetrads){
	foreach my $tetrad (sort keys %{$tetrads{$plate}}){
		foreach my $strain (sort keys %{$tetrads{$plate}{$tetrad}}){
			print ASSESS "$strain\t$plate\t$tetrad\t$strain_info{$strain}[3]\t$assess{$strain}[0]\t$assess{$strain}[1]\t$assess{$strain}[3]";
			foreach my $chromosome (sort keys %positions){
				unless (${$assess{$strain}[2]}{$chromosome}){${$assess{$strain}[2]}{$chromosome}="-";}
				print ASSESS "\t${$assess{$strain}[2]}{$chromosome}";
			}
			print ASSESS "\n";
		}
	}
}
close ASSESS;

##################################################################################
############################# Assess Marker Metrics ##############################
##################################################################################

#Assess each marker for the number of strains with parent 1 allele call, parent 2
#allele call and the error calls 0 and 3. Store counts in hash with chromosome
#, position and allele type (0,1,2,3 etc) as keys to count of that allele call for 
#that marker across all strains in the dataset using the confidence method (not 
#the perfect match counting method)
#e.g.$positions{chr1}{1001}{1}=5, is 5 counts of allele 1 at chr1 pos 1001
&assess_markers(\%positions, \%markers, "$outdir/initial_marker_assessment.tsv", $label);

##################################################################################
############################## Delete Bad Markers ################################
##################################################################################

#Delete any bi-alleic markers called less than 50 times or where the ratio of P1 
#allele calls to P2 alleles calls is less than 4/5 or more than 5/4 and any 
#mono-allelic calls occurring in more than 60% of all strains
&delete_bad_markers(\%markers, \%positions, $outdir);

##################################################################################
############### Include Dead Spores in their Correct Tetrads #####################
##################################################################################

&add_dead_spores(\%markers, \%strain_info, \%tetrads, "$outdir/added_dead_strains.tsv", $dead_file);

##################################################################################
############################# Call Missing Markers ###############################
##################################################################################

#So that all markers are called in all strains, call missing markers as zero in all
#strains 
&set_missing_to_zero(\%markers, \%strain_info, \%positions);

##################################################################################
########################### Reinitialize Bad Strains #############################
##################################################################################

#Set all allele calls from high (>10%) heterozygosity strains to 0, do the same 
#thing to strains that appear to be the single incorrectly assigned strain in a 
#tetrad
&reinitialize_bad_strains(\%markers, \%assess, \%tetrads, \%strain_info, $outdir);

##################################################################################
############################# Tetrad Inference ###################################
##################################################################################

#Use tetrad information to infer missing allele calls.
#If 3 other strains have a P1 or P2 allele call can always infer the missing fourth.
#If 2 other strains have both P1 calls or both P2 calls then can infer the 
#missing values for the two missing strains. Note that this assumes the called
#values are correct and assigns no confidence to the inferred values
&infer_by_tetrad_segregation(\%markers, \%tetrads, \%positions, \%strain_info);

##################################################################################
##################### Calculate Genetic Map and Distances ########################
##################################################################################

&calculate_recombination_map(\%markers, "$outdir/genetic_map.tsv", "$outdir/genetic_distances.tsv");

##################################################################################
############ Calculate Relationship Between Physical Distance and cM #############
##################################################################################

my $gen_coeff=&map_physical_to_genetic_distance("$outdir/gen.r", "$outdir/genetic_distances.tsv", "$outdir/gen_coeff.tsv", $r_path);
#my $gen_coeff=3.645562e-06;
#my $tet_coeff=3.645562e-06;


##################################################################################
####################### Calculate Tetrad Map and Distances #######################
##################################################################################

#Tetrad map and calculating of tetrad coefficient not used as very slow and 
#coefficient should be identical to that calculated from single spore recombination
#frequencies

#&calculate_tetrad_map(\%tetrads, \%markers, \%positions, "$outdir/tetrad_map.tsv", "$outdir/tetrad_distances.tsv");

##################################################################################
#### Calculate Relationship Between Physical Distance and Tetrad Transitions #####
##################################################################################

#my $tet_coeff=&map_physical_to_tetrad_distance("$outdir/tet.r", "$outdir/tetrad_distances.tsv", "$outdir/tet_coeff.tsv", $r_path);
my $tet_coeff=$gen_coeff;

##################################################################################
############ Infer Missing Data Using Tetrad Transition Probabilities ############
##################################################################################

&infer_by_tetrad_map(\%tetrads, \%markers, \%positions, \%strain_info,  $tet_coeff);

##################################################################################
############## Infer Missing Data Using Recombination Probabilities ##############
##################################################################################

&infer_by_genetic_map(\%markers, \%positions, \%strain_info, $gen_coeff);

##################################################################################
##################### Reassess Markers After Inference Steps ####################
##################################################################################

&assess_markers(\%positions, \%markers, "$outdir/final_marker_assessment.tsv", $label);

##################################################################################
##################### Reassess Genetic Map and Distances #########################
##################################################################################

&calculate_recombination_map(\%markers, "$outdir/genetic_map_post_inference.tsv", "$outdir/genetic_distances_post_inference.tsv");

##################################################################################
############################# Output Haplotypes ##################################
##################################################################################

&p1_p2_support(\%markers, "$outdir/final_p1_p2_support.tsv");


#Output final combined haplotype file
open (OUT, ">$outdir/hap.tsv") or die "Could not open inferred haplotype file $outdir/hap.tsv\n$!\n";
#print OUT "#Strain\tChromosome\tPosition\tBarcode\tPlate\tTetrad\tSpore\tScore_P1\tScore_P2\tSupporting_Reads*\tAllele\tPerfect_P1_reads\tPerfect_P2_reads\tTotal_perfect_reads\tAllele_call_2\tYSEQ\tCross\n";
#print OUT "#*If inferred then method and Phred confidence (TI-tetrad 2:2 inference, TP-tetrad recombination pattern inference, RC-spore linkage inference). No confidence value for TI.\n";
print OUT "strain\tchromosome\tposition\tbarcode\tplate\ttetrad\tspore\tscore_p1\tscore_p2\tsupporting_reads\tallele\tperfect_p1_reads\tperfect_p2_reads\ttotal_perfect_reads\tallele_call_2\tyseq\tcross\n";
foreach my $strain (sort keys %markers){
	foreach my $chromosome (sort keys %{$markers{$strain}}){
		foreach my $position (sort {$a<=>$b} keys %{$markers{$strain}{$chromosome}}){
		my $line=join ("\t", @{$markers{$strain}{$chromosome}{$position}});
		print OUT "$strain\t$chromosome\t$position\t$line\n";
		}
	}
}
	
close OUT;




##################################################################################
##################################################################################
############################### Subroutines ######################################
##################################################################################
##################################################################################

#Adds dead spores from lims YPG strain csv output file that are members of tetrads
#sequenced by radtag back into the strain-marker, strain-info and tetrad hashes.
#Sets strain-marker hash only at the strain level, later code will add marker info
#Strain-info set with plate and tetrad info from csv file to values as normal except
#barcode and yseq are set to "DEAD"
#Writes added strains to file, with strain, ycr, plate tetrad and spore for the added
#inviable strain
#Arguments:Pointers to strain-marker, strain-info and tetrad hashes, output file and
#csv formatted input file from LIMS, with strain, ycr, plate and tetrad as columns
#1,7,8 and 12.
sub add_dead_spores{

	my ($markers, $strain_info, $tetrads, $outfile, $deadfile)=@_;

	open (DEAD, "<$deadfile") or die "Could not open file $deadfile describing non viable strains:\n$!\n";
	open (OUT, ">$outfile") or die "Could not write file $outfile listing non viable strains:\n$!\n";
	print OUT "#Strain\tCross\tPlate\tTetrad\tSpore\n";

	my @spores=qw(A B C D);
	while (<DEAD>){
		chomp;
		my @elements=split /,/;

		#Find dead spores
		unless ($elements[-1] eq "BLANK"){next;}
		my ($strain, $ycr, $plate, $tetrad)=@elements[0,6,7,11];
		
		#Check that strain is associated with existing tetrad
		unless ($$tetrads{$plate}{$tetrad}){next;}

		#Find other members of tetrad
		my @other_strains=keys %{$$tetrads{$plate}{$tetrad}};

		#Identify spores already set
		my %seen_spores=();
		foreach my $other_strain (@other_strains){
			$seen_spores{$$strain_info{$other_strain}[3]}=1;
		}
		
		#Assign first unused spore value to this strain
		#Add strain to hashes and output file descfibing added dead strains
		foreach my $spore (@spores){
			if ($seen_spores{$spore}){next;}
			else {
				print OUT "$strain\t$ycr\t$plate\t$tetrad\t$spore\n";
				@{$$strain_info{$strain}}=("DEAD", $plate, $tetrad, $spore, "DEAD", $ycr);
				$$tetrads{$plate}{$tetrad}{$strain}=1;
				%{$$markers{$strain}}=();
				last;
			}

		}

	}
	close OUT;
	close DEAD;

}

#Counts the number of markers is each strain called as P1 or P2 and outputs to file for all strains
#Arguments:Pointer to strain-marker hash and output filename
sub p1_p2_support{

	my ($markers, $outfile)=@_;

	#Hash holding counts of each observed allele call for all strains:strain as key to array
	#of total markers, total markers with P1 or P2 call
	my %counts=();

	foreach my $strain (keys %$markers){
		@{$counts{$strain}}=(0,0);
		foreach my $chromosome (keys %{$$markers{$strain}}){
			foreach my $position (keys %{$$markers{$strain}{$chromosome}}){
				++$counts{$strain}[0];	#Count all markers
				if ($$markers{$strain}{$chromosome}{$position}[7]==1 or $$markers{$strain}{$chromosome}{$position}[7]==2){++$counts{$strain}[1]}	#Count if P1 or P2
			}
		}
	}

	#Write numbers to output file
	open (OUT, ">$outfile") or die "Could not write to file $outfile:\n$!\n";

	print OUT "#Strain\ttotal_strain_markers\tP1_P2_markers\n";
	foreach my $strain (sort keys %counts){
		print OUT "$strain\t$counts{$strain}[0]\t$counts{$strain}[1]\n";
	}
	close OUT;
}


#Ensures that all marker is position hash are called in all strains by setting missing markers in
#marker-position hash to zero for all count, support and allele array positions.Note that plate,
#spore, yseq and ycr values in marker hash are set correctly from strain-info hash
#Arguments: pointers to strain-marker, strain info and marker-position hashes.

sub set_missing_to_zero{

	my ($markers, $strain_info, $positions)=@_; 

	#Check all strains in strain-marker hash for all marker positions in positions hash
	#if absent initilaize with 0 values
	foreach my $chromosome (keys %$positions){
		foreach my $position (keys %{$$positions{$chromosome}}){

			foreach my $strain (keys %$strain_info){
				unless ($$markers{$strain}{$chromosome}{$position}){
					@{$$markers{$strain}{$chromosome}{$position}}=(@{$$strain_info{$strain}}[0,1,2,3],0,0,0,0,0,0,0,0,$$strain_info{$strain}[4], $$strain_info{$strain}[5]);
				}
			}
		}
	}
}

#Calculates inferred values across a set of contiguous markers in a tetrad using the 
#allele calls that do exist in the tetrad along with converting the physical distance
#to genetic distance
#Arguments:Coefficient to convert physical to genetic distance, pointers to reduced tetrad,
#additions hashes and current chromosome. Additions hash has chromosome, marker and spore as
#key to array where element 1 is probability of allele1 across all posible tetrad recombination
#patterns, elements 2 is prob of allele 2, element 3 is allele call and element 4 is phred value
#for the called allele. Reduced tetrads hash is chromosome and position as key to array of allele
#calls (1,2, or 0 for all non 1 or 2 calls) ordered by strain in the current tetrad.
sub ditype_inference{

	my ($coeff, $tetrad, $additions, $chromosome)=@_;

	#Set the six possible tetrad states
	my %states=(
	"0"=>[1,1,2,2],
	"1"=>[1,2,1,2],
	"2"=>[1,2,2,1],
	"3"=>[2,2,1,1],
	"4"=>[2,1,2,1],
	"5"=>[2,1,1,2],
	);


	#Record emission "probabilities" for each state- for each spore in tetrad
	#will either be 0 or 1 depending on current state. i.e. for state 0 (pattern 1,1,2,2)
	#the porbability of spore 1 emission being 1 is 1 and the probability of it being 2 is zero
	my %e_probs=();


	#Set the emission probabilities for each state
	#State, spore and parental allele as key to 1 or 0 probabilities
	foreach my $state (keys %states){
		foreach my $index (0..3){
			if ($states{$state}[$index]==1){
				$e_probs{$state}{$index}{1}=1;
				$e_probs{$state}{$index}{2}=0;
			} else {
				$e_probs{$state}{$index}{1}=0;
				$e_probs{$state}{$index}{2}=1;

			}
		}
	}


	#Calculate the number of recombination events separating each pair of states- should be 0 (PD), 1 (TT) or 2 (NPD)
	#is equal to the number of crossover spres divided by 2
	#First state and second state as key to number of crossovers
	my %transition_steps=();

	foreach my $first_state (keys %states){
		foreach my $second_state (keys %states){
			my $differences=0;
			foreach (0..3){
				unless ($states{$first_state}[$_]==$states{$second_state}[$_]){++$differences;}
			}

			$transition_steps{$first_state}{$second_state}=$differences/2;
		}
	}	

	
	#Initialize variables
	my %current_transitions;	#Probability of PD, NPD, TT for current physical interval. Number of cos (0,2,1) as key to probability
	my $current_pos=0;		#Position being considered
	my $old_pos=0;			#Last position
	my @positions=sort {$a<=>$b} keys %$tetrad;	#ordered positions from contiguous region in which to carry out inference
	my %paths=();			#State path across the positions. Key is the set of ordered states for each called position (e.g. 01231) as key to path probability

	#Do not attempt to calculate all of the paths across more than 8 markers-takes too long/too much memory
	if (@positions>8){return;}

	foreach my $index (0..$#positions){

		#Get current position and distance from last position, convert into genetic distance
		$current_pos=$positions[$index];
		my $dist=$current_pos-$old_pos;
		my $co_count=$dist*$coeff;	#Convert physical distance into Morgans using coefficient calculated earlier

		#Calculate the transition probabilities given the current genetic distances for the three types of transition: 
		#parental ditype, tetratype and non parental ditype
	
		#Tetra-type
		$current_transitions{1}=2/3*(1-exp(-3*$co_count));

		#There are 4 tetratype patterns so the probability of any single pattern is 1/4 the total tetratype probability
		$current_transitions{1}=$current_transitions{1}/4;

		#PD (only 1 PD pattern)
		$current_transitions{0}=1/6*(1+2*exp(-3*$co_count)+3*exp(-2*$co_count));

		#NPD (only 1 NPD pattern)
		$current_transitions{2}=1/6*(1+2*exp(-3*$co_count)-3*exp(-2*$co_count));

		#First position- ignore distance from last position and transition proabilities, only calculate initial emission probabilities for paths
		#Starting with every possible state at this position (removing those that have emission proability of zero)
		if ($index==0){

			#Inialise 6 paths, each starting with one of the possible tetrad marker states, calculate emission probabilities
			foreach my $start (keys %states){
				my $probability=1;
				foreach my $spore (0..3){
					if ($$tetrad{$current_pos}[$spore]==0){next}
					$probability*=$e_probs{$start}{$spore}{$$tetrad{$current_pos}[$spore]};
				}
				unless ($probability){next;}	#Ignore initial states that are not possible- do not build paths from these
				$paths{$start}=$probability;
			}
		}

		#Internal positions use distance from last marker to calculate transition porbabilities for each possible state at the new position and combine
		#with emission probability for that state at that position remove paths getting zero probability (impossible emission pattern)
		else {
			#Extend each exiting path with each of the 6 possible states at this position
			foreach my $path (keys %paths){
				if (length($path)<$index){delete $paths{$path}; next;}			#Unextended paths shorter than last position- remove any such (is this necessary?)
				my $old_probability=$paths{$path};	#probability of path to be extended
				my @nodes=split ('', $path);		#path sequence of states
				my $last_node=$nodes[-1];		#last state in path

				#Foreach possible new state....
				foreach my $next_state (keys %states){

					#calculate probability of transition from old to new state and multiply with original path probability
					my $probability=$old_probability * $current_transitions{$transition_steps{$last_node}{$next_state}};

					#Update probability with emission probability for current state (will be 0 or 1)
					foreach my $spore (0..3){
						if ($$tetrad{$current_pos}[$spore]==0){next}
						$probability*=$e_probs{$next_state}{$spore}{$$tetrad{$current_pos}[$spore]};
					}
					unless ($probability){next;}	#Ignore new paths that are not possible- do not build further paths from these
					$paths{$path.$next_state}=$probability;	#Register new path and its probability
				}
				#Remove path from hash after adding the 6 possible new extended versions of the path
				delete $paths{$path};
			}
		}
		$old_pos=$current_pos;	#Update position
	}


	#Remove any short or zero proability paths
	foreach my $path (keys %paths){
		if (length($path)<=$#positions){delete $paths{$path}; next;}
		unless ($paths{$path}){delete $paths{$path}; next;}
	}


	#For each 0 position in original tetrad calculate probabilities of P1 and P2 across all remaining possible paths, add the probability of each path to the
	#allele it supports that this position:chromosome position and spore as key to array with index 1 = support for allele 1 and index 2 = support for allele 2
	foreach my $index (0..$#positions){
		my $position=$positions[$index];
		foreach my $spore (0..3){
			if ($$tetrad{$position}[$spore]==0){
				foreach my $path (keys %paths){
					my @path_array=split ('', $path);
					$$additions{$chromosome}{$position}{$spore}[$states{$path_array[$index]}[$spore]]+=$paths{$path};
				}
			}
		}
	}



	foreach my $position (sort {$a<=>$b} keys %{$$additions{$chromosome}}){
		foreach my $spore (sort {$a<=>$b} keys %{$$additions{$chromosome}{$position}}){

			#if no support for allele 1 or no support for allele 2 set the probabilities to zero
			unless ($$additions{$chromosome}{$position}{$spore}[1]){$$additions{$chromosome}{$position}{$spore}[1]=0;}
			unless ($$additions{$chromosome}{$position}{$spore}[2]){$$additions{$chromosome}{$position}{$spore}[2]=0;}

			#If no support for allele 1 then probability of allele 2 is infinitely small
			if ($$additions{$chromosome}{$position}{$spore}[1]==0 and $$additions{$chromosome}{$position}{$spore}[2]>0){
				$$additions{$chromosome}{$position}{$spore}[3]=2;
				$$additions{$chromosome}{$position}{$spore}[4]="INF";
			}

			#If no support for allele 2 then probability of allele 1 is infinitely small
			elsif ($$additions{$chromosome}{$position}{$spore}[2]==0 and $$additions{$chromosome}{$position}{$spore}[1]>0){
				$$additions{$chromosome}{$position}{$spore}[3]=1;
				$$additions{$chromosome}{$position}{$spore}[4]="INF";
			}

			#If support for both the probability of each is its support/sum of support for both alleles. Calculate and convert to phred (-10Log10 P value)
			elsif ($$additions{$chromosome}{$position}{$spore}[1]>$$additions{$chromosome}{$position}{$spore}[2]){
				$$additions{$chromosome}{$position}{$spore}[3]=1;
				$$additions{$chromosome}{$position}{$spore}[4]=-10*log($$additions{$chromosome}{$position}{$spore}[2]/($$additions{$chromosome}{$position}{$spore}[1] +$$additions{$chromosome}{$position}{$spore}[2]))/log(10);
				$$additions{$chromosome}{$position}{$spore}[4]=int($$additions{$chromosome}{$position}{$spore}[4]+0.5);
			}
			elsif	($$additions{$chromosome}{$position}{$spore}[2]>$$additions{$chromosome}{$position}{$spore}[1]){
				$$additions{$chromosome}{$position}{$spore}[3]=2;
				$$additions{$chromosome}{$position}{$spore}[4]=-10*log($$additions{$chromosome}{$position}{$spore}[1]/($$additions{$chromosome}{$position}{$spore}[1] +$$additions{$chromosome}{$position}{$spore}[2]))/log(10);
				$$additions{$chromosome}{$position}{$spore}[4]=int($$additions{$chromosome}{$position}{$spore}[4]+0.5);
			}
		}
	}
}


#Use known markers and physical distances between markers to be inferred and known 
#markers to infer the missing values. Relies on global estimate of how physical 
#distance converts to genetic distance 
#Arguments:Pointers to strain-marker,position hashes and strain-info hashes, coefficient to convert 
#physical distance (bp) to genetic distance (M)
sub infer_by_genetic_map{
	
	my ($markers, $all_positions, $strain_info, $coeff)=@_;

	#Infer markers one strain and one chromosome at a time
	foreach my $strain (keys %$markers){
		foreach my $chromosome (keys %$all_positions){

			my %rec_infer=();	#Holds inferred markers to be added to strain-marker hash

			#Identify all missing calls for this strain on this chromosome
			my @missing=();
			my @temp=(); #Get plate tetrad etc info about this strain
			unless ($$markers{$strain}{$chromosome}){next;}	#If no markers on this chromosome then cannot infer missing values
			foreach my $position (sort {$a<=>$b} keys %{$$all_positions{$chromosome}}){
				#Attempt to infer all non P1 or P2 allele calls
				unless ($$markers{$strain}{$chromosome}{$position} and ($$markers{$strain}{$chromosome}{$position}[7]==1 or $$markers{$strain}{$chromosome}{$position}[7]==2)){push (@missing, $position);}
			}
			foreach my $position (@missing){
				#For each missing call find the first leftward and rightward calls that are not inferred by tetrad recombination
				my @left=();
				my @right=();
				#Examine all positions on chromosome with allele call, get the P1 or P2 calls closest to left and right side of marker to be inferred
				foreach my $second_position (sort {$a<=>$b} keys %{$$markers{$strain}{$chromosome}}){
					if  ($$markers{$strain}{$chromosome}{$second_position}[6]=~/^TP/){next;}	#Ignore inferred markers for doing further inference
					unless  ($$markers{$strain}{$chromosome}{$second_position}[7]==1 or $$markers{$strain}{$chromosome}{$second_position}[7]==2){next;}
					
					#Record each position as leftward until position is greater than marker being inferred. Take its allele call and its distance from the marker to infer.
					elsif ($second_position<$position){@left=($$markers{$strain}{$chromosome}{$second_position}[7], $position-$second_position);}
					#Take first position past marker to be inferred as rightward marker.  Take its allele call and its distance from the marker to infer.
					elsif ($second_position>$position){@right=($$markers{$strain}{$chromosome}{$second_position}[7], $second_position-$position); last;}
				}
	
				#Initialise confidence and inferred allele call
				my $phred=0;
				my $inferred=0;

				#Both left and right flanking markers identified
				if (@left and @right){

					#Calculate probabilty of left and right crossovers using Haldane's function with genetic distance inferred from physical distance using coefficient
					my $p_left=(1-exp(-2*$coeff*$left[1]))/2;
					my $p_right=(1-exp(-2*$coeff*$right[1]))/2;
					
					#Flanks in agreement with each other: only double co or zero co possible, call as flank allele
					if ($left[0]==$right[0]){
						$inferred=$left[0];
						$phred=-10*log($p_left*$p_right/($p_left*$p_right+(1-$p_left)*(1-$p_right)))/log(10)	#only double co or zero co possible, error if double
					}
				
					#Flanks not in agreement, marker closer to left flank therefore will be called as left flank allele
					elsif ($left[1]<$right[1]){
						$inferred=$left[0];
						$phred=-10*log($p_left/(1-$p_left*$p_right-(1-$p_left)*(1-$p_right)))/log(10);	#only single leftward or rightward crossovers possible, error if left
					}
					#Flanks not in agreement, marker closer to right flank therefore will be called as right flank allele
					elsif ($right[1]<$left[1]){
						$inferred=$right[0];
						$phred=-10*log($p_right/(1-$p_left*$p_right-(1-$p_left)*(1-$p_right)))/log(10);	#only single leftward or rightward crossovers possible, error if right

					}
				}

				#Marker on only left flank, call as left flank
				elsif (@left){
					$inferred=$left[0];
					$phred=-10*log((1-exp(-2*$coeff*$left[1]))/2)/log(10);	#probability incorrect is frequency of left crossover
				}				
				
				#Marker on only right flank, call as right flank
				elsif (@right){
					$inferred=$right[0];
					$phred=-10*log((1-exp(-2*$coeff*$right[1]))/2)/log(10);	#probability incorrect is frequency of right crossover
				}	

				#No flanking  markers
				else {next;}

				#Round phred score
				$phred=int($phred+0.5);
				if ($phred<20){next;}	#Only report inference if confidence better than phred 19.
				@{$rec_infer{$strain}{$chromosome}{$position}}=(@{$$strain_info{$strain}}[0,1,2,3],0,0,"RC_${phred}", $inferred,0,0,0,0, @{$$strain_info{$strain}}[4,5]);
					

			}

			#Apply inferred values for this chromosome to strain-marker hash 	
			foreach my $position (keys %{$rec_infer{$strain}{$chromosome}}){
 				@{$$markers{$strain}{$chromosome}{$position}}=@{$rec_infer{$strain}{$chromosome}{$position}};
			}
		}
	}
}

#Use known markers across each tetrad and physical distances between markers 
#to be inferred and known markers to infer the missing values. Relies on global 
#estimate of how physical distance converts to genetic distance and hence different
#bimarker tetrad patterns (PD, TT, NPD) 
#Arguments:Pointers to tetrad, strain-marker,position hashes and strain-info hashes, 
#coefficient to convert physical distance (bp) to genetic distance (M)
sub infer_by_tetrad_map{

	my ($tetrads, $markers, $all_positions, $strain_info, $coeff)=@_;

	foreach my $plate (keys %$tetrads){
		foreach my $tetrad (keys %{$$tetrads{$plate}}){

			#Get array of strain names for current tetrad
			my @strains=sort keys %{$$tetrads{$plate}{$tetrad}};	

			#Ignore bad tetrads
			unless (@strains==4 and $$tetrads{$plate}{$tetrad}{$strains[0]}){next;}

			#Move through the genome one chromosome at a time 
			#finding places where there are not 4 values called and attempting to infer the missing values 
			my $last=0;	#Holds previous position considered
			my $no_capture=1;
			foreach my $chromosome (sort keys %$all_positions){

				my %additions=();#holds inferred values to be included in strain-marker hash. Format strain chromosome marker as keys to array
						#where index 1 is P1 probability, index 2 is p2 probability, index 3 is allele call and index 4 is phred support for allele call
				my %current=();	#Holds array of allele calls across the members of the tetrad with position as key, note all non 1 or 2 calls are converted to zeros in array (will be inferred)
						#Acumulates consecutive missing value positions and any single flanking fully called positions. Inferred values calculated from the available data for each such
						#contiguous region	
				foreach my $position (sort {$a<=>$b} keys %{$$all_positions{$chromosome}}){
					foreach my $strain (@strains){
						if ($markers{$strain}{$chromosome}{$position} and ($$markers{$strain}{$chromosome}{$position}[7]==1 or $$markers{$strain}{$chromosome}{$position}[7]==2)){
							push (@{$current{$position}},$$markers{$strain}{$chromosome}{$position}[7]); 
						} else{
							push (@{$current{$position}},0); 
	
						}
					
					}
				
					#Full data at this position (no zero calls) and last position also fully called (no_capture flag still set), then delete previous position 
					if (!(join ('',@{$current{$position}})=~/0/) & $no_capture){
						delete $current{$last};
					}
	
					#Incomplete region at end of chromosome, i.e. current position has missing values and is last position on the chromosome
					#Use data in %current chromosome, which has all of the positions with missing data back to first fully called position
					#to call the missing values.
					elsif (((join ('',@{$current{$position}}))=~/0/) & ($position==(sort {$a<=>$b} keys %{$$all_positions{$chromosome}})[-1])){
						&ditype_inference($coeff, \%current, \%additions, $chromosome);
						foreach my $back_pos (keys %current){
							unless ($back_pos==$position){delete $current{$back_pos};} 
						}
						my $no_capture=1;
			
					}
						
		
					#Starting or continuing an incomplete region
					elsif ((join ('',@{$current{$position}}))=~/0/){
						$no_capture=0;
					}
	
					#Full data after incomplete region infer the missing values held in %current.
					#Reinitialise %current after completion
					else {
						&ditype_inference($coeff, \%current, \%additions, $chromosome);
						foreach my $back_pos (keys %current){
							unless ($back_pos==$position){delete $current{$back_pos};} #Keep final position as may be needed to infer next set of missing values
						}
						my $no_capture=1;
											
					}
					$last=$position;	#Update previous position
				
				}

				#Add the inferred values to the strain-position hash if the confidence of the inferred values is at least Phred 20
				foreach my $position (sort keys	%{$additions{$chromosome}}){		
					foreach my $spore (0..3){
						if ($additions{$chromosome}{$position}{$spore}[3]){
							my $inferred=$additions{$chromosome}{$position}{$spore}[3];
							my $phred=$additions{$chromosome}{$position}{$spore}[4];
	
							my $strain=$strains[$spore];
	
							#Set the fields for current strain-marker combo to inferrence 0 counts for P1 and P2 support and "TP" and Phred for total reads then inferred value as haplotype
							#WILL ONLY APPLY IF THIS STRAIN HAS A MISSING (non P1 or P2) VALUE AT THIS MARKER
							unless ($phred eq "INF" or $phred>=20){next;}
							@{$$markers{$strain}{$chromosome}{$position}}=(@{$$strain_info{$strain}}[0,1,2,3],0,0,"TP_${phred}", $inferred,0,0,0,0, @{$$strain_info{$strain}}[4,5]);
	
						}
					}					
				}
			}
		}
	}
}



#Uses table of tetrad transition frequencies (NPD, PD, TT) by marker distance to calculate 
#relationship between physical distance and average number of crossovers. Returns 
#coefficient to translate physical distance (bp) into Morgans
#Arguments: filename for r commands, genetic distance table (6 columns, physical distance between
#markers and frequency of 0,1,2,3,4 positions crossed over (0=PD, 2=TT, 4=NPD, 1 and 3 should
#not occur), coefficient output file from R, path to R executable
sub map_physical_to_tetrad_distance{

	my ($rfile, $infile, $outfile, $r_path)=@_;

	#Create R batch file to estimate conversion of physical distance to M
	#using Haldane's tetrad mapping function
	open (R_BATCH, ">$rfile") or die "Could not open R batch file $rfile:\n$!\n";


	print R_BATCH "tetrads<-read.table(\"$infile\");\n";
	print R_BATCH "temp<-tetrads[tetrads\$V4<300000,];\n";
	print R_BATCH "tt_test=nls(temp\$V7~2/3*(1-exp(-3*a*temp\$V4)), start=list(a=1/300000), data=temp);\n";
	print R_BATCH "pd_test=nls(temp\$V5~1/6*(1+2*exp(-3*a*temp\$V4)+3*exp(-2*a*temp\$V4)), start=list(a=1/300000), data=temp);\n";
	print R_BATCH "t=mean(c(as.numeric(coef(tt_test)), (as.numeric(coef(pd_test)))));\n";
	print R_BATCH "write (t, \"$outfile\");\n";

	close R_BATCH;

	#Execute R script and write out conversion coefficient, extract coefficient from
	#this file and return it
	system "$r_path CMD BATCH $rfile";

	open (COEFF, "<$outfile") or die "Could not open coefficient file $outfile\n$!\n";
	my $coeff=<COEFF>;
	chomp $coeff;
	close COEFF;
	return $coeff;
}

#Uses table of rf values by marker distance to calculate relationship between physical distance and
#average number of crossovers. Returns coefficient to translate physical distance (bp) into Morgans
#Arguments: filename for r commands to be outputted, genetic distance table (2 tab delimited columns,
# physical distance between markers and rf), coefficient output filename from R, path to R executable
sub map_physical_to_genetic_distance{


	my ($rfile, $infile, $outfile, $r_path)=@_;

	#Create R batch file to estimate conversion of physical distance to M
	#using Haldane's mapping function
	open (R_BATCH, ">$rfile") or die "Could not open R batch file $rfile:\n$!\n";

	print R_BATCH "temp<-read.table(\"$infile\");\n";
	print R_BATCH "recoms<-temp[temp\$V1<300000,];\n";
	print R_BATCH "test=nls(recoms\$V2~(1-exp(-2*a*recoms\$V1))/2, start=list(a=1/300000), data=recoms);\n";
	print R_BATCH "write (t, \"$outfile\");\n";

	close R_BATCH;

	#Execute R script and write out conversion coefficient, extract coefficient from
	#this file and return it
	system "$r_path CMD BATCH $rfile";

	open (COEFF, "<$outfile") or die "Could not open coefficient file $outfile\n$!\n";
	my $coeff=<COEFF>;
	chomp $coeff;
	close COEFF;

	return $coeff;
}

#Takes a set of allele calls for markers from tetrads, considering tetrads one
#by one for each pair of marker positions and makes a strain specific pairwise hash 
#indicating whether the pair show a non-parental ditype, tetratype or parental ditype
#recombination pattern. Then updates counts of these patterns for the marker pair 
#across all tetrads. Once all tetrads processed calculates proportion of each pattern
#seen in all tetrads for that marker pair. Outputs file with distance between markers
#as first column, then frequency of difference counts 0-4 (representing PD, errror,
#TT, error and NPD patterns)
#Arguments: Pointers to tetrad hash, strain-marker hash, position hash and output filenames
#for tetrad pattern map (not currently used) and distance-pattern file.
sub calculate_tetrad_map  {

	my ($tetrads, $markers, $all_positions, $outfile1, $outfile2)=@_;

	#Holds counts of frequency of PD, TT and NPD patterns for each pair of markers on each chromosome
	#across all strains. Chromosome, first marker and second marker as keys to array where index gives
	#number of differences (0=PD, 2=TT, 4=NPD, 1 and 3 are errors) holding counts of number of times
	#that number of differences/pattern seen. Index 5 gives total counts,	
	my %difference_matrix=();

	foreach my $plate (keys %$tetrads){
		foreach my $tetrad (keys %{$$tetrads{$plate}}){
			unless ($$tetrads{$plate}{$tetrad}{(keys %{$tetrads{$plate}{$tetrad}})[0]}){next;}	#Ignore bad tetrads for calculation


			#Get array of strain names for current tetrad
			my @strains=keys %{$tetrads{$plate}{$tetrad}};


			#Move through the genome finding places where there are 4 values called that are all P1 or P2  
			foreach my $chromosome (sort keys %$all_positions){
				my %current=();			
				foreach my $first_position (sort {$a<=>$b} keys %{$$all_positions{$chromosome}}){
					my @first_array=();	#Array of values for current marker position for current tetrad
					my $bad_flag=0;
					foreach my $strain (sort @strains){
						if ($$markers{$strain}{$chromosome}{$first_position} and ($$markers{$strain}{$chromosome}{$first_position}[7]==1 or $$markers{$strain}{$chromosome}{$first_position}[7]==2)){
							push (@first_array,$$markers{$strain}{$chromosome}{$first_position}[7]);
						} else{
							$bad_flag=1;
	
						}
					
					}
	
					#Ignore current position if P1 or P2 not called in all 4 strains
					if ($bad_flag){next;}
					unless (@first_array==4){next;}


					#Compare tetrad pattern to all LATER tetrads on the same chromosome	
					foreach my $second_position (sort {$a<=>$b} keys %{$$all_positions{$chromosome}}){
						if ($first_position>=$second_position){next;}
						my @second_array=();	#Array of values for this second marker position
						my $bad_flag2=0;
						foreach my $strain (sort @strains){
							if ($$markers{$strain}{$chromosome}{$second_position} and ($$markers{$strain}{$chromosome}{$second_position}[7]==1 or $$markers{$strain}{$chromosome}{$second_position}[7]==2)){
								push (@second_array,$$markers{$strain}{$chromosome}{$second_position}[7]);
							} else{
								$bad_flag2=1;
							}
					
						}
						#Ignore position if P1 or P2 not called in all 4 strains
						if ($bad_flag2){next;}
						unless (@second_array==4){next;}
	
	
						my $differences=0;	#Number of spores in tetrad with recombination 0=PD, 2=TT, 4=NPD
						#If current marker pair not considered for previous tetrad initialize it now (first position, second position as keys to counts of number of difference (as index in array)
						unless ($difference_matrix{$chromosome}{$first_position}{$second_position}){@{$difference_matrix{$chromosome}{$first_position}{$second_position}}=(0,0,0,0,0,0);}
				
						foreach my $index (0..3){
							unless ($first_array[$index]==$second_array[$index]){++$differences;}
						}
	
						#Update counts of the tetrad bilocus pattern (PD, TT, NPD etc) for this pair of loci and update total count (element[5])
						++$difference_matrix{$chromosome}{$first_position}{$second_position}[$differences];
						++$difference_matrix{$chromosome}{$first_position}{$second_position}[5];
	
	
					}
				}
			}
		}
	}

	#For each pair of positions on each chromosome output the proportion of PD, TT, NPD and error (1 or 3 differences) calls for each marker pair on each chromosome
	#In form of distance between markers (first column) and then proportion of each number of differences as columns 2-6
	open (TETRADS, ">$outfile2") or die "Could not open tetrad recombination frequncey file $outfile2:\n$!\n";

	foreach my $chromosome (keys %difference_matrix){
		foreach my $first_pos (keys %{$difference_matrix{$chromosome}}){
			foreach my $second_pos (keys %{$difference_matrix{$chromosome}{$first_pos}}){
				unless ($difference_matrix{$chromosome}{$first_pos}{$second_pos}[5]>10){next;}	#Only calculate frequencies when minimum of 10 comparisons
				#Calculate frequency of each pattern type
				my $zero_rat=$difference_matrix{$chromosome}{$first_pos}{$second_pos}[0]/$difference_matrix{$chromosome}{$first_pos}{$second_pos}[5];
				my $first_rat=$difference_matrix{$chromosome}{$first_pos}{$second_pos}[1]/$difference_matrix{$chromosome}{$first_pos}{$second_pos}[5];
				my $second_rat=$difference_matrix{$chromosome}{$first_pos}{$second_pos}[2]/$difference_matrix{$chromosome}{$first_pos}{$second_pos}[5];
				my $third_rat=$difference_matrix{$chromosome}{$first_pos}{$second_pos}[3]/$difference_matrix{$chromosome}{$first_pos}{$second_pos}[5];
				my $fourth_rat=$difference_matrix{$chromosome}{$first_pos}{$second_pos}[4]/$difference_matrix{$chromosome}{$first_pos}{$second_pos}[5];
	
	
				print TETRADS "$chromosome\t";
				print TETRADS "$first_pos\t";
				print TETRADS "$second_pos\t";
	
				print TETRADS$second_pos-$first_pos."\t".(int($zero_rat*100+0.5))/100;
				print TETRADS"\t".(int($first_rat*100+0.5))/100;
				print TETRADS "\t".(int($second_rat*100+0.5))/100;
				print TETRADS "\t".(int($third_rat*100+0.5))/100;
				print TETRADS "\t".((int($fourth_rat*100+0.5))/100)."\n";
			}
		}
	}

	close TETRADS;
}

#Takes a set of allele calls for markers from a single strain and makes a strain 
#specific pairwise hash indicating whether the pair have alleles from the same 
#parent ,or different (recombinant). Then updates count of each marker pair seen 
#and whether recombinant in  a second hash that includes all strains and can be used 
#to calcuate pairwise rf values across the strains. Recmbination hash has first chrom, 
#first pos, second chrom, second pos as key to array of total comparisons, total 
#recombinant comparisons. This hash is then outputted to a recombination table
#with markers as rows and columns (first row/column=chromosome, second row/column=
#position) values in table are rf values. Also outputs a file with 2 columns giving
#the physical distance between each pair of markers on the same chromosome and the rf
#betweem them.
#Arguments: Pointers to strain-marker hash, rf table file, physical distance to rf file
sub calculate_recombination_map{

	my ($current_markers, $outfile1, $outfile2)=@_;

	my %all_markers=();

	#Convert allele calls for markers in current strain into a pairwise hash with first chromosome and position, second chromosome and position as keys to an array
	#where the first value is 1 (indicating this pairwise comparison exists in the current strain) and then either 0 if same parental allele at both loci or 1 if alleles at
	#the two loci come from different parents recombinant)
	#For each marker
	foreach my $strain (keys %$current_markers){
		my %current_map=();	#Hash for current strain
		foreach my $first_chromosome (keys %{$$current_markers{$strain}}){
			foreach my $first_position (keys %{$$current_markers{$strain}{$first_chromosome}}){
			 	my $first_allele=$$current_markers{$strain}{$first_chromosome}{$first_position}[7]; 	#Identify parental allele at first position
				unless ($first_allele==1 or $first_allele==2){next;}	#Only P1 and P2 alleles are useful for calculating recombination

				#Compare to all other markers
				foreach my $second_chromosome (keys %{$$current_markers{$strain}}){
					foreach my $second_position (keys %{$$current_markers{$strain}{$second_chromosome}}){
						my $second_allele=$$current_markers{$strain}{$second_chromosome}{$second_position}[7]; #Identify parental allele at second position
						unless ($second_allele==1 or $second_allele==2){next;}	#Only P1 and P2 alleles are useful for calculating recombination
						++$current_map{$first_chromosome}{$first_position}{$second_chromosome}{$second_position}[0];	#Register current pair

						
						#If not recombinant then set second array position to 0, otherwise 1 
						if ($first_allele==$second_allele){
							$current_map{$first_chromosome}{$first_position}{$second_chromosome}{$second_position}[1]=0;
						}

						else {$current_map{$first_chromosome}{$first_position}{$second_chromosome}{$second_position}[1]=1;}
					}
				}
			}
		}

		#Update pairwise comparisons for all strains using info from current strain
		foreach my $first_chromosome (keys %current_map){
			foreach my $first_position (keys %{$current_map{$first_chromosome}}){
				foreach my $second_chromosome (keys %{$current_map{$first_chromosome}{$first_position}}){
					foreach my $second_position (keys %{$current_map{$first_chromosome}{$first_position}{$second_chromosome}}){
						
						#Update total number of times pair seen
						$all_markers{$first_chromosome}{$first_position}{$second_chromosome}{$second_position}[0]+=$current_map{$first_chromosome}{$first_position}{$second_chromosome}{$second_position}[0];
						#Update number of time pair seen with non parental allele configuration
						$all_markers{$first_chromosome}{$first_position}{$second_chromosome}{$second_position}[1]+=$current_map{$first_chromosome}{$first_position}{$second_chromosome}{$second_position}[1];
					}
				}
			}
		}
	}

	#Output rf table to file and output table of physical distance-rf pairs for calculating conversion of physical distance to genetic distance (same chromosome marker comparisons only
	open (RECOM, ">$outfile1") or die "Could not open recombination map file $outfile1:\n$!\n";
	open (RECOM_DIST, ">$outfile2") or die "Could not open distance-rf file $outfile2:\n$!\n";

	
	#Print column headers:chromosome in row 1, position in row 2
	print RECOM "\t";
	foreach my $first_chromosome (sort keys %all_markers){
		foreach my $first_position (sort {$a<=>$b} keys %{$all_markers{$first_chromosome}}){
			print RECOM "\t$first_chromosome";
		}
	}
	print RECOM "\n";


	print RECOM "\t";
	foreach my $first_chromosome (sort keys %all_markers){
		foreach my $first_position (sort {$a<=>$b} keys %{$all_markers{$first_chromosome}}){
			print RECOM "\t$first_position";
		}
	}
	print RECOM "\n";

	#For each maker print a row, with chromosome, then position of marker then rf to each other marker by column
	foreach my $first_chromosome (sort keys %all_markers){
		foreach my $first_position (sort {$a<=>$b} keys %{$all_markers{$first_chromosome}}){
			print RECOM "$first_chromosome\t$first_position";
			foreach my $second_chromosome (sort keys %{$all_markers{$first_chromosome}{$first_position}}){
				foreach my $second_position (sort {$a<=>$b} keys %{$all_markers{$first_chromosome}{$first_position}{$second_chromosome}}){
					my $rat=int(($all_markers{$first_chromosome}{$first_position}{$second_chromosome}{$second_position}[1]/$all_markers{$first_chromosome}{$first_position}{$second_chromosome}{$second_position}[0]) * 1000 +0.5)/1000;

					print RECOM "\t",$rat;

					#If same chromosome output distance-rf pair
					if ($first_chromosome eq $second_chromosome and $second_position>$first_position){
						print RECOM_DIST abs($second_position-$first_position)."\t$rat\n";
					}
				}
			}
			print RECOM "\n";
		}
	}

	close RECOM_DIST;
	close RECOM;

}



#Uses tetrad segregation pattern (2:2 for any marker- ignoring gene conversions)
#to infer missing data. If 3/4 strains have P1/P2 calls then the 4th can always
#be inferred. If 2/4 have allele calls the other 2 can be inferred 1/2 time.
#Makes no attempt to measure confidence of prediction. Inferred row has 0
#support values, has the P1/P2 confidence allele call set and has "TI" (tetrad inference)
#in the reads supporting confidence call column.

#Additionally identifies occurences of 3 P1 or 3 P2 allele calls at marker positions
#within tetrads and sets all values for these over-represented calls to 0. Note 
#might be better to reset the lowest confidence of these calls.
#Arguments:Pointers to strain-marker, tetrad, position and strain info hashes.
sub infer_by_tetrad_segregation {

	my ($markers, $tetrads, $positions, $strain_info)=@_;

	#Examine each tetrad, one at a time
	foreach my $plate (sort keys %$tetrads){
		foreach my $tetrad (sort keys %{$$tetrads{$plate}}){

			#Get all strains in this tetrad
			my @all_strains=keys %{$$tetrads{$plate}{$tetrad}};
			unless (@all_strains){next;}
			#Ignore tetrads that show a lot of 3:1 and 4:0 segregation behaviour i.e. set to 0 in tetrad hash
			if ($$tetrads{$plate}{$tetrad}{$all_strains[0]}){

				#Check segregation pattern across all positions for this tetrad, missing values will be set as 0
				foreach my $chromosome (sort keys %$positions){
					foreach my $position (sort {$a<=>$b} keys %{$$positions{$chromosome}}){
	

						#Get allelle calls for these strains at the current position, if marker absent then set as zero
						#Initialise occurances of allele 0,1,2 to zero
						my %parent_calls=();
						$parent_calls{0}=0;
						$parent_calls{1}=0;
						$parent_calls{2}=0;
						my @calls=();
						my @sample_line=();
						foreach my $strain (@all_strains){
							if (defined $$markers{$strain}{$chromosome} and defined $$markers{$strain}{$chromosome}{$position} and defined ($$markers{$strain}{$chromosome}{$position})){
								++$parent_calls{$$markers{$strain}{$chromosome}{$position}[7]};	#Increment count of current allele
								push (@calls, $$markers{$strain}{$chromosome}{$position}[7]);	#Push current allele onto tetrad array
							}
							else {++$parent_calls{0}; push (@calls, 0);}
						
						}

						#If two calls of allele 1, the other two strains must be allele 2
						if ($parent_calls{1}==2){
							foreach my $index (0..$#all_strains){
								#Set non-allele 1 calls to allele 2 in strain-marker hash
								unless ($calls[$index]==1){
									my $strain=$all_strains[$index];
									@{$$markers{$strain}{$chromosome}{$position}}=(@{$$strain_info{$all_strains[$index]}}[0,1,2,3],(0,0,"TI",2,0,0,0,0),$$strain_info{$strain}[4],$$strain_info{$strain}[5]);
								}
							}
						}

						#If two calls of allele 2, the other two strains must be allele 1
						elsif ($parent_calls{2}==2){
							foreach my $index (0..$#all_strains){
								#Set non-allele 2 calls to allele 1 in strain-marker hash
								unless ($calls[$index]==2){
									my $strain=$all_strains[$index];
									@{$$markers{$strain}{$chromosome}{$position}}=(@{$$strain_info{$all_strains[$index]}}[0,1,2,3],(0,0,"TI",1,0,0,0,0),$$strain_info{$strain}[4],$$strain_info{$strain}[5]);
								
								}
							}
						}


						#If 1 allele called more than twice in tetrad, reinitialise the too-frequent calls to zero and infer by other methods	
						elsif ($parent_calls{1}>2){
							foreach my $index (0..$#all_strains){
								if ($calls[$index]==1){
									my $strain=$all_strains[$index];
									@{$$markers{$strain}{$chromosome}{$position}}=(@{$$strain_info{$all_strains[$index]}}[0,1,2,3],(0,0,0,0,0,0,0,0),$$strain_info{$strain}[4],$$strain_info{$strain}[5]);
								
								}
							}

						}
						#If 2 allele called more than twice in tetrad, reinitialise the too-frequent calls to zero and infer by other methods	
						elsif ($parent_calls{2}>2){
							foreach my $index (0..$#all_strains){
								if ($calls[$index]==2){
									my $strain=$all_strains[$index];
									@{$$markers{$strain}{$chromosome}{$position}}=(@{$$strain_info{$all_strains[$index]}}[0,1,2,3],(0,0,0,0,0,0,0,0),$$strain_info{$strain}[4],$$strain_info{$strain}[5]);
								
								}
							}

						}
					}
				}
			}
		}
	}
}


#Takes strains that have high heterozygosity (>10%) or which appear to be the single incorrect strain in a 
#tetrad and reinitialise all marker calls to allele 0 in the strain-marker hash. Record in "reinitialized strains.tsv"
#file. Rescue tetrads if only one problem strain and good segregation among rest of markers
#Arguments"Pointers to strain-marker hash, assemment hash, tetrad hash and strain-info hash, output directory 
sub reinitialize_bad_strains{

	my ($markers, $assess, $tetrads, $strain_info, $outdir)=@_;

	my %resets=();	#Hold strains to be reset
	open (RESETS, ">$outdir/reinitialized strains.tsv") or die "Could not write to renitialized strain file $outdir/reinitialized strains.tsv:\n$!\n";

	#Flag highly heterozygous strains for resetting and check tetrads with low coverage or highly heterozygous strains for
	#good marker segregation among the other memebers of the tetrad, if true reset tetrad to "good"
	foreach my $strain (sort keys %$markers){
		if ($$strain_info{$strain}[4] eq "DEAD"){next;}	#No assessment for dead strains???

		#Only interested in highly heterozygous or low coverage strains
		unless ($$assess{$strain}[1]>10 or $$assess{$strain}[0]<1000){next;}

		#Flag highly heterozygous strains for reset
		if ($$assess{$strain}[1]>10){
			print RESETS "Resetting\t$strain\theterozygosity too high\n";
			$resets{$strain}=1;
		}
		
		#Reset tetrad to good if 3:0 segregation rare in the rest (non-heterozygous or non-low coverage) of the members of tetrad
		if ($$assess{$strain}[3] eq "NA"){next;}	#No assessemnt of other members of the tetrad
		elsif ($$assess{$strain}[3]<=10){
			my ($plate, $tetrad)=@{$$strain_info{$strain}}[1,2];
			foreach my $all_strain (keys  %{$$tetrads{$plate}{$tetrad}}){
				$$tetrads{$plate}{$tetrad}{$all_strain}=1;
			}
		}			
	}

	#Identify strains incorrectly assigned to tetrad (i.e. appears 1/4 wrong), identify for reset and rescue tetrad
	foreach my $plate (keys %$tetrads){
		foreach my $tetrad (keys %{$$tetrads{$plate}}){
			my @selected=();
			my $not_selected='';
			foreach my $strain (keys %{$$tetrads{$plate}{$tetrad}}){
				if ($$tetrads{$plate}{$tetrad}{$strain}){next;}	#Ignore good tetrads
				#if (!($$assess{$strain}[3] eq "NA") and $$assess{$strain}[3]>10) {push (@selected, $strain);}
				#Place strains that have NA or >10% het assessment into array for this tetrad. Any low het strain will be set to "not_selected"
				if ($$assess{$strain}[3] eq "NA" or $$assess{$strain}[3]>10) {push (@selected, $strain);}
				else {$not_selected=$strain;}
			}

			#Rescue tetrad if only 3 strains show high het/NA and 4th shows low het			
			if (@selected==3 and $not_selected){
				$resets{$not_selected}=1;
				print RESETS "Resetting\t$not_selected\tnot member of $plate $tetrad\n";
				foreach my $strain (keys %{$$tetrads{$plate}{$tetrad}}){
					$$tetrads{$plate}{$tetrad}{$strain}=1;
				}
			}
		}
	}

	#Reset to 0 all markers in strains with too high heterozygosity or where placed
	#in wrong tetrad
	foreach my $strain (sort keys %resets){
		foreach my $chromosome (keys %{$$markers{$strain}}){
			foreach my $position (keys %{$$markers{$strain}{$chromosome}}){
				@{$$markers{$strain}{$chromosome}{$position}}[4,5,6,7,8,9,10,11]=(0,0,0,0,0,0,0,0);
			}
		}
	}
	close RESETS;
}
			

#Removes markers with suspicious segregation patterns from the positions and strain-marker hashes
#Bi-allelic markers are removed if the allele1/allele2 ratio is smaller than 0.8 or greater than 1.25
#Mono-alleleic markers are removed if allele1 or allele2 calls occur in >60% of all strains (this is not
# very effective filter). Strains with fewer than 50 total P1 and P2 calls cannot be assessed for
#correct segregation if bi-allelic and are also removed
# Outputs a file "deleted_markers.tsv" to the output directory listing deleted 
#markers in format chromosome, position, bi-allelic/mono-allelic and reason for removal.
#Arguments: Pointers to strain-marker and position hashes and output directory for deleted markers file
sub delete_bad_markers{

	my ($markers, $positions, $outdir)=@_;

	#Identify markers that pass the quality tests
	#Marker must be called as 1,2,4,5,6
	my $total_strains=scalar (keys %$markers);
	my %good_markers=();	#Holds chromosome and position of markers passing the quality tests

	my $ratio=0.8;	#Minimum ratio when number strains with less frequent parental allele call is divided by number of strains with more frequent
		#identifies bi-alleic markers with strange segregation ratio

	my $mono_limit=0.6;	#Maximum proportion of strains that can have a mono-allelic call (should be 0.5 if called in all strains where expected)
	my $minimum_calls=50;

	#Create file to hold deleted markers failing quality control
	open (DEL, ">$outdir/deleted_markers.tsv") or die "Could not write to deleted markers file $outdir/ideleted_markers.tsv:\n$!\n";

	print DEL "#Chromosome\tPosition\tMarker_type\tReason_for_Deletion\n";
 
	#Assess each bi-allelic marker
	foreach my $chromosome (sort keys %$positions){
		foreach my $position (sort {$a<=>$b} keys %{$$positions{$chromosome}}){
			
			unless ($$positions{$chromosome}{$position}{"type"} eq "BA"){next;}	#Only consider bi-allelic markers	
			#Record any marker in good marker hash as long as meets ratio requirements
			my ($smaller_count, $greater_count)=sort {$a<=>$b} ($$positions{$chromosome}{$position}{1}, $positions{$chromosome}{$position}{2});
			my $coverage=$smaller_count + $greater_count;
			unless ($coverage>$minimum_calls){print DEL "$chromosome\t$position\tBi-Allelic\ttoo few P1 and P2  calls\n"; next;}		#Fails test
			unless ($smaller_count){print DEL "$chromosome\t$position\tBi-Allelic\tratio zero\n"; next;}					#Fails test
			if (($smaller_count/$greater_count) < $ratio){print DEL "$chromosome\t$position\tBi-Allelic\tratio outside 0.8-1.25\n";}	#Fails test,
			else {$good_markers{$chromosome}{$position}=1;}											#Passes test, record
		}
	}

	#Assess each mono-allelic marker
	foreach my $chromosome (keys %$positions){
		foreach my $position (sort {$a<=>$b} keys %{$$positions{$chromosome}}){
			
			unless ($$positions{$chromosome}{$position}{"type"} eq "MA"){next;}	#Only consider mono-allelic markers	
			#Record any marker in good marker hash as long as not called as allele 1 (or 2) in more than 60% of all strains
			my ($greater_count)=(sort {$a<=>$b} ($$positions{$chromosome}{$position}{1}, $positions{$chromosome}{$position}{2}))[1];
			if (($greater_count/$total_strains) > $mono_limit ){print DEL "$chromosome\t$position\tMono-Allelic\tcalled in >60% of strains\n";}	#Fails test
			else {$good_markers{$chromosome}{$position}=1;}												#Passes test, record
		}
	}

	close DEL;

	
	#Delete badly any marker not recorded as being good from strain-marker and position hashes
	foreach my $strain (sort keys %$markers){
		foreach my $chromosome (sort keys %{$$markers{$strain}}){
			foreach my $position (sort {$a<=>$b} keys %{$$markers{$strain}{$chromosome}}){
				unless ($good_markers{$chromosome}{$position}){

					#Remove from both hashes
					delete $$markers{$strain}{$chromosome}{$position};
					delete $$positions{$chromosome}{$position};

					
					#unless ($$positions{$chromosome}){delete $positions{$chromosome};}
					#unless ($$markers{$strain}{$chromosome}) {delete $$markers{$strain}{$chromosome};}
					#unless ($$markers{$strain}) {delete $$markers{$strain};}
				}
			}
		}
	}
}


#Counts occurences of each allele call (0,1,2,3 etc) for each marker across all strains
#Populates positions hash with chromosome, position and allele as keys to counts of that
#allele. Outputs marker assessment table with chromosome and position columns, then columns
#of counts for each allele call. Tab delimited "marker_assessment.tsv" file. 
#Arguments: Pointers to positions and strain-marker hashes, output assessment file, label
#to identify cross-iteration
sub assess_markers {
	my ($positions, $markers, $outfile, $label)=@_;

	#Zero out all counts of markers
	foreach my $chromosome (sort keys %$positions){
		foreach my $position (sort keys %{$$positions{$chromosome}}){
			foreach (0..6){
				$$positions{$chromosome}{$position}{$_}=0;
			}
		}
	}

	#Move through strain-marker hash and count occurances of all allele types for all marker positions
	foreach my $strain (sort keys %$markers){
		foreach my $chromosome (sort keys %{$$markers{$strain}}){
			foreach my $position (sort {$a<=>$b} keys %{$$markers{$strain}{$chromosome}}){
				
				my $code=$$markers{$strain}{$chromosome}{$position}[7];	#Current allele
				$$positions{$chromosome}{$position}{$code}+=1;		#Update allele count for this position
			}
		}
	}

	#Count total number of strains
	my $total_strains=keys %markers;


	#Output marker assessment file:each marker with total number of strains in which it occurs, and counts of each code-0 to 6 for that
	#marker
	open (OUT, ">$outfile") or die "Could not open output file $outfile\n";

	#print OUT "#Chromosome\tposition\ttotal_strains\tcode0_strains\tcode1_strains\tcode2_strains\tcode3_strains\tcode4_strains\tcode5_strains\tcode6_strains\n";
	print OUT "chromosome\tposition\ttotal_strains\tcode0_strains\tcode1_strains\tcode2_strains\tcode3_strains\tcode4_strains\tcode5_strains\tcode6_strains\cross_id\n";
	foreach my $chromosome (sort keys %$positions){
		foreach my $position (sort {$a<=>$b} keys %{$$positions{$chromosome}}){
			print OUT "$chromosome\t$position\t$total_strains";
			for (my $i=0; $i<7; ++$i){
				unless ($$positions{$chromosome}{$position}{$i}){$$positions{$chromosome}{$position}{$i}=0;}	#Create counts of 0 for codes never observed for this strain
				print OUT "\t$$positions{$chromosome}{$position}{$i}";
			}
			print OUT "$label\n";
		}
	}
	close OUT;
}



#For each strain, see if the other 3 members of the tetrad have any 3:0 segregation this will identify any single
#misclassified strain as it will be the only strain member where this frequency is close to 0, the others will have higher values.
#Report this value for all strains as 4th element of assessment hash array for that strain. Also assesses if any member of a 
#tetrad has >10% 3:0 segregation among the rest of the tetrad. If so sets the tetrad as "bad", in tetrad hash all strains 
#have associated value 0, rather than 1
#Arguments:Pointers to assessment, tetrad and strain-marker hashes
sub assess_tetrads{
	
	my ($assess, $tetrads, $markers, $positions)=@_;

	foreach my $plate (keys %$tetrads){
		foreach my $tetrad (keys %{$$tetrads{$plate}}){

			#Analyse all of the spores in the current tetrad

			#Get array of strain names for current tetrad
			my @strains=keys %{$$tetrads{$plate}{$tetrad}};

			#For each strain assess frequency of 3:0 allele segregation  across all markers
			#among the other strains in the tetrad
			#only count markers where P1 and P2 allele calls were made in the 3 other strains 
			foreach my $strain (@strains){
				my @calls=(0,0);	#Count of aberrant (3:1 or 4:0 segregation) and 2:2 segregation across all analyzable markers
				foreach my $chromosome (sort keys %$positions){
					foreach my $position (sort {$a<=>$b} keys %{$$positions{$chromosome}}){
						my @temp=(0,0,0,0,0,0,0);	#Counts of each haplotype within the tetrad (haplotype is index)
	

						#Examine alleles for other member sof the tetrad
						foreach my $other(@strains){
							unless (defined($$markers{$other}{$chromosome}{$position})){next;}	#Marker not present
							my $current_allele=$$markers{$other}{$chromosome}{$position}[7];
							if (!($strain eq $other)){++$temp[$current_allele];}	#Do not compare current strain to itself
						}

						if ($temp[1]+$temp[2]==3){++$calls[1];}	#Each of the other strains in the tetrad have a P1 or P2 call
						if ($temp[1]==3 or $temp[2]==3){++$calls[0];}	#All of the other strains have EITHER a P1 or a P2 call for this marker (aberrant)

					}	
				}
			
				#If no tetrads counted (none having 3 calls of P1+P2 among the 3 other members of the tetrad) then cannot calculate proportion having aberrant segregation
				if ($calls[1]==0){$$assess{$strain}[3]="NA";}

				#Record frequency of 3 P1 or 3 P2 calls among other memebers of tetrad for this strain. If too high- note that the tetrad has aberrant segregation behaviour
				else{
					$$assess{$strain}[3]=int(($calls[0]/$calls[1]*100)+0.5); 

					#If more than 10% of markers typed show 3:0 segregation among the other strains for any chosen strain
					#mark this tetrad by setting 0 in place of 1 for value in tetrad hash for all strains in this tetrad
					if (int(($calls[0]/$calls[1]*100)+0.5)>10){
						$$tetrads{$plate}{$tetrad}{$_}=0 foreach @strains;	
					}
				}


			}
			
		}
	}
}

#Assess each strain for quality metrics:strain as key to array:
#Total reads, % informative markers called as allele3 (likely mixed sample), 
#hash holding ploidy for each chromosome (not 2 micron or mitochondria),
#percentage 3:0 calls of P1 or P2  among the other 3 members of the same tetrad,
#where all calls P1 or P2 (single low value when rest of tetrad high 
#indicates single spore incorrectly assigned to tetrad)
#Arguments:Assessment hash pointer, tetrad hash pointer, strain-marker hash pointer, marker position hash pointer
sub assess_strains{

	my ($assess, $tetrads, $markers, $positions)=@_;

	#Assess strains 1 at a time	
	foreach my $strain (sort keys %$markers){

		my %chromosome_count=();	#Hash to hold first read counts and then ploidy for this strain across all chromosome
						#Chromosome as key to read counts or ploidy for that chromosome

		my %calls=();		#Temp hash to hold counts of each allele call type for this strain (0,1,2,3 etc) for calculating "heterozygosity":allele call as key to count of occurance
		my $total_calls=0;	#Total positions with non 0 allele call for this strain

		#Process marker information for this strain 1 marker at a time
		foreach my $chromosome (sort keys %{$$markers{$strain}}){
			foreach my $position (sort {$a<=>$b} keys %{$markers{$strain}{$chromosome}}){
				++$calls{$$markers{$strain}{$chromosome}{$position}[7]};		#Increment counts of current allele call for this strain
				$$assess{$strain}[0]+=$markers{$strain}{$chromosome}{$position}[6];	#Increment read counts for this strain
				$chromosome_count{$chromosome}+=$$markers{$strain}{$chromosome}{$position}[6];		#Increment read counts for this chromosome


				unless ($markers{$strain}{$chromosome}{$position}[7]==0){++$total_calls;}	#Increment called markers for this strain (unless marker 0)
				
			}
		}

		#Record by chromosome read counts for this strain in assessment hash
		$$assess{$strain}[2]=\%chromosome_count;
	
		#Calculate % of markers called (haplotypes 1, 2, 3 (not 0)) that are 'heterozygous' (haplotype 3) and record
		if ($total_calls and $calls{3}){$$assess{$strain}[1]=int(($calls{3}/$total_calls*100)+0.5);}	
		else {$$assess{$strain}[1]=0;}	#No non 0 calls, means no haplotype 3 calls
	}

	#Convert the read counts by chromosome across all strains to ploidy estimates
	&convert_counts_to_ploidy($assess);

	#Examine marker behaviour in each tetrad (for each strain freq of 3:0 segregation among the other members of its tetrad)
	#result will be item 4 in hash array
	&assess_tetrads($assess, $tetrads, $markers, $positions);


}


#Converts proportions of all reads in a strain found on each chromosome into a ploidy assessment
#using the first quartile proportion of reads aligned for that chromosome across all strains
#Argument: Assessment hash pointer where strain is key to an array, the 1st element of which is the
#total number of aligned reads for that strain and the 3rd element of which is
#a pointer to a hash with chromosome name as key to number of reads aligned to that chromosome
sub convert_counts_to_ploidy{
	
	my ($assess)=@_;


	#Calculate for each strain, the proportion of all reads that come from each chromosome, do not include 2micron or mitochondrial genomes
	#(removed earlier)	
	my %ratios=();
	foreach my $strain(sort keys %$assess){
		foreach my $chromosome (sort keys %{$$assess{$strain}[2]}){

			#Proportion of all reads for this strain that come from the current chromosome
			my $ratio=${$$assess{$strain}[2]}{$chromosome}/$assess{$strain}[0];
			$ratios{$chromosome}{$strain}=$ratio;
		}
	}

	#For each chromosome, look across all strains and identify the proportion of all genomic reads that come from that chromosome
	# Calculate the 1st quartile of these values. This can then be used to classify each strain as dismic/monosomic for that chromosome. 
	#Assume at most 50% of strains disomic, hence first quartile will be in the distribution of ratio values for monosomic strains.
	#Hold quartile values in hash %quarts with chromosome as key
	my %quarts=();
	foreach my $chromosome (keys %ratios){
		my @proportion=();	#For current chromosome, all of the read proportions coming from this chromosome across all strains
		foreach my $strain (keys %{$ratios{$chromosome}}){
			push (@proportion, $ratios{$chromosome}{$strain});
		}
		@proportion=sort {$a<=>$b} @proportion;	#order the array of values
		my $quart=int(scalar(@proportion)/4);	#Index of first quartile value in the sorted array
		$quart=$proportion[$quart-1];	#-1 as 0 based indexing, get the 1st quartile proportion
		$quarts{$chromosome}=$quart;	#Set current chromosome to first quartile value
	}
	
	#Normalize  each strain chromosome combo ratio by the first quartile value for that chromosome
	#to convert to ploidy value estimate
	foreach my $strain (keys %$assess){
		foreach my $chromosome (keys %{$$assess{$strain}[2]}){
			${$$assess{$strain}[2]}{$chromosome}=int(($ratios{$chromosome}{$strain}/$quarts{$chromosome})*10+0.5)/10;
		}
	}

}


#Combine all of the strain-yseq haplotype calls for all strains in the cross directories
#Convert haplotype call file into a hash with YPG number (strain), chromosome and position as keys
#to array of tag, plate, tetrad, spore, p1 confidence, p2 confidence, total informative reads, allele call,
#,P1 perfect match reads, P2 perfect match reads, perfect match allele call, yseq number and cross number (YCR)
#Also populates tetrad, position and strain info hashes:
#
#Tetrads: plate, tetrad and strain as keys to value 1 if tetrad has consistent marker segregation (good) or 0 if not 
#(set to 1 here)
#
#Positions:Chromosome, position and word "type" as keys to BA (bi-allelic- parent specific reads expected 
#from both parents) or MA (mono-allelic, restriction site present in only one parent)
#
#Strain info:strain as key to array of tag, plate, tetrad,spore, YSEQ and YCR (cross)
#
#Note assumes that a single strain is not sequenced more than once- if it is then the highest numbered
#YSEQ version will be the only one kept.
#Arguments:Hash pointers to hold information by strain+marker, by tetrad,by marker and
#by strain and pointer to array of yseq folder names containing the raw haplotype files 
sub combine_yseqs{

	my ($markers, $tetrads, $positions, $strain_info, $yseq_folders)=@_;


	my %no_duplicates=();	#Strain and yseq as keys to value 1, use to detect strains with more than 1 yseq (error)
	foreach my $yseq_folder (@$yseq_folders){

		#Get current ycr
		$yseq_folder=~/\/(YCR\d+)\//;
		my $ycr=$1;

		#Get raw haplotype file for current YSEQ
		open (IN, "<$yseq_folder/raw_haplotypes.txt") or die "Could not open haplotype file $yseq_folder/raw_haplotypes.txt:\n$!\n";

		#Get YSEQ number
		my $yseq=basename $yseq_folder;

		#Process haplotype file from this YSEQ
		while (<IN>){
			chomp;
			if (/^#/){next;}		#Ignore comment line

			#Extract columns from raw haplotype files
			my @elements=split /\t/;


			#Ignore non- progeny lines
			unless ($elements[1]){next;}	#Ignore tags with no associated YPG or YO number
			$elements[1]=uc($elements[1]);		#Uppercase names

			if ($elements[1]=~/^YO/){next;}		#Ignore YO names- only interested in progeny here
			if ($elements[1]=~/^YPG490/ or $elements[1]=~/^YPG512/){next;}		#Ignore YO names- only interested in progeny here

			#Ignore sequences for 2-micron or mitochondria
			if ($elements[5] eq "chr2m" or $elements[5] eq "chrmt"){next;}

			#Ignore lines with 0 support (caused by reads with ambiguous alignment)
			if ($elements[7] eq "NA"){unless ($elements[8]>0){next;}}
			elsif ($elements[8] eq "NA"){unless ($elements[7]>0){next;}}
			else {unless ($elements[7]>0 or $elements[8]>0){next;}}

			#Record tetrad info:plate, tetrad and strain as keys to value 1
			$$tetrads{$elements[2]}{$elements[3]}{$elements[1]}=1;

			#Set current strain-marker info: strain, chromosome and position as keys to array (see above)
			@{$markers{$elements[1]}{$elements[5]}{$elements[6]}}=(@elements[0,2,3,4,7,8,9,10,11,12,13,14], $yseq, $ycr);

			#Keep track of strain information:strain as key to tag, plate, tetrad, spore, YSEQ and YCR
			@{$$strain_info{$elements[1]}}=($elements[0], $elements[2], $elements[3], $elements[4], $yseq, $ycr);

			#Record this position (chromosome and position as keys to bi-allelic, mono-allelic call)
			$$positions{$elements[5]}{$elements[6]}{"type"}="BA";
			if ($elements[7] eq "NA" or $elements[8] eq "NA"){$$positions{$elements[5]}{$elements[6]}{"type"}="MA";}

			#Test that all strains have been sequenced only once. Always take highest numbered yseq if duplicates.
			$yseq=~/YSEQ(\d+)/;
			$no_duplicates{$elements[1]}{$1}=1;
			if ((keys %{$no_duplicates{$elements[1]}})>1){
				print "Duplicate sequencing of strain $elements[1]:\n";
				until ((keys %{$no_duplicates{$elements[1]}})==1){
					my $lowest=(sort keys %{$no_duplicates{$elements[1]}})[0];
					delete $no_duplicates{$elements[1]}{$lowest};
					print "\t\t\t\t\tDeleting YSEQ${lowest}\n";
				}
				
			}
		}
		close IN;
	}
}


